
Completion of IGN data, road example
Source:vignettes/web_only/completion_of_ign_data_road_example.Rmd
completion_of_ign_data_road_example.Rmd
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(happign) # retrieve pnr boundary
#> Warning in CPL_gdal_init(): GDAL Message 1: Unable to find driver DODS to
#> unload from GDAL_SKIP environment variable.
#> Warning in CPL_gdal_init(): GDAL Message 1: Unable to find driver DODS to
#> unload from GDAL_SKIP environment variable.
#> Please make sure you have an internet connection.
#> Use happign::get_last_news() to display latest geoservice news.
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tmap);tmap_mode("view") # interactive map
#> tmap mode set to interactive viewing
library(osmextract)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
#> Check the package website, https://docs.ropensci.org/osmextract/, for more details.
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, unionContext
To compare the completeness of IGN data, OpenStreetMap data will be used. They will be retrieved with the package osmextract.
The study area is the territory regional natural park of the Massif des Bauges located in the Northern Alps, in the Bauges massif. It has an area of about 900 km2, spread over the departments of Savoie and Haute-Savoie. It is made up of 67 communes.
Extract PNR boundary
To extract IGN data with happign, it is necessary to
have a geographical input. I chose the commune of Annecy located in the
PNR of the Massif des Bauges. happign contains the
cog_2022 dataset which allows to find the insee code of
this commune to retrieve boundary with
get_apicarto_commune.
code_insee <- cog_2022[cog_2022$LIBELLE == "Annecy", 1][1] #use cog_2022 happign dataset to find insee code
annecy <- get_apicarto_commune(code_insee) # get annecy boundary
pnr <- get_wfs(annecy,
apikey = "environnement",
layer_name = "PROTECTEDAREAS.PNR:pnr") # search for pnr containinng annecy
qtm(pnr) #quick checkExtract roads data
## From ign bdcarto with happign
all_road_ign <- get_wfs(pnr,
apikey = "cartovecto",
layer_name = "BDCARTO_BDD_WLD_WGS84G:troncon_route")
road_ign <- st_intersection(all_road_ign, pnr)
## from osm data with osmextract
all_road_osm <- oe_get("Savoie",
provider = "openstreetmap_fr",
stringsAsFactors = FALSE,
force_download = TRUE,
force_vectortranslate = TRUE,
quiet = TRUE)
road_osm <- st_intersection(all_road_osm, pnr) |>
filter(highway %in% c("primary", "primary_link", "secondary","secondary_link",
"tertiary", "tertiary_link", "trunk", "trunk_link",
"residential", "cycleway", "living_street", "unclassified",
"motorway", "motorway_link", "pedestrian", "steps",
"track","service", NA))Results
The data set provided by osmextract contains 4 times more linear than the BDCARTO of the ign. Be careful however, the osmextract data requires more verifications than the IGN data : 25% of the OSM dataset are unqualified roads (highway == NA).
# who has the longest ?
sum(st_length(road_ign))
#> 1683201 [m]
sum(st_length(road_osm))
#> 6773427 [m]
# plot results
map <- tm_shape(pnr)+
tm_borders()+
tm_shape(road_osm)+
tm_lines(col = "firebrick", lwd = 2)+
tm_shape(road_ign)+
tm_lines(col = "#95C019", lwd = 2)+
tm_add_legend(type = "fill",
labels = c("IGN Roads", "OSM Roads"),
col = c("#95C019", "firebrick"))
map